## pipeline script for MethyPattern align in to genome feature
## author: Yaping Liu  lyping1986@gmail.com
my $mem="4g";
my $data_dir = "/export/uec-gs1/laird/shared/research/yaping/code/GNOME_seq/MethyPatternFeatureWalker/HCT116_hg19/MethyPatternResult/promoter/";
my $feature_dir = "/export/uec-gs1/laird/shared/research/yaping/code/GNOME_seq/MethyPatternFeatureWalker/HCT116_hg19/location_data/";
my $plot_dir = "/export/uec-gs1/laird/shared/research/yaping/code/GNOME_seq/MethyPatternFeatureWalker/HCT116_hg19/plot/promoter/";
my $genome_assembly = "hg19";
my $bissnp_dir = "/home/uec-00/yapingli/software/BisSNP/BisSNP-0.68.4.jar";
my $R_plot_script = "/home/uec-00/yapingli/code/mytools/R/MethyPatternFeaturePlot.R";
#my $bam_file = "/home/uec-00/yapingli/data/Terry/fullBam/157-NS_merge_D09AYACXX_D0ECCACXX_KEL656A79.hg19_rCRSchrm.fa.mdups.reheader.bam";
my $bam_files_dir = "/export/uec-gs1/laird/shared/production/ga/flowcells/C0J2KACXX/run1/results/C0J2KACXX/";

my $reference_seq;
my $dbsnp;
if($genome_assembly eq "hg19"){
	$reference_seq = "/export/uec-gs1/laird/shared/research/yaping/data/genome_data/genome/hg19/hg19_rCRSchrm.fa";
	$dbsnp = "/export/uec-gs1/laird/shared/research/yaping/data/genome_data/dbSNP/dbsnp_135.hg19.sort.vcf";
}
else{
	$reference_seq = "/export/uec-gs1/laird/shared/research/yaping/data/genome_data/genome/hg18/hg18_unmasked.plusContam.fa";
	$dbsnp = "/export/uec-gs1/laird/shared/research/yaping/data/genome_data/dbSNP/dbsnp_135.hg18.sort.vcf";
}


my $data_matrix_scale = 2000;
my $confidance_level = 20;
my $feature_name = 1;
my $plot_step = 20;
my $plot_scale = 1000;
my $plot_axis_step = 500;
my $smooth = 1; ## 0: no smooth; 1: smooth

my $count = 1;
my $use_age = "USAGE: perl MethyPatternFeaturePipeline.pl data_dir feature_dir plot_dir";
#if($#ARGV <3){
#	print "$use_age\n";
#	exit(1);
#}

opendir(DH1,$bam_files_dir);
	foreach my $sub_dir(readdir(DH1)){
		next if($sub_dir =~ /\./ or $sub_dir =~ /\.\./);
		my $bam_file = $bam_files_dir."$sub_dir/"."ResultCount_".$sub_dir.".hg19_rCRSchrm.fa.mdups.bam";
		my $sample="FL_".$sub_dir;
		#print "$sub_dir\n";
		opendir(DH,$feature_dir);
		foreach my $feature_file(readdir(DH)){
		if($feature_file =~ /(\S+)\.minusUpstream0\.plusDownstream0(\S+sort.bed)$/){
			my $prefix = $1;
			my $suffix = $2;
			my $region_file = $prefix.".minusUpstream3000\.plusDownstream3000".$suffix;
			my $header = "#PBS -q laird\n";
			$header .= "#PBS -l walltime=168:00:00,mem=$mem,nodes=1:ppn=4\n";
			#$header .= "#PBS -l walltime=168:00:00,mem=$mem\n";
			$header .= "#PBS -N $feature_file-$sample\n";
			$header .= "#PBS -j oe\n";
			$header .= "cd \$PBS_O_WORKDIR\n";
			my $java_cmd = "java -Xmx$mem -jar $bissnp_dir -T MethyPatternFeature ";
			$java_cmd .= "-R $reference_seq -I $bam_file -D $dbsnp -stand_call_conf $confidance_level -stand_emit_conf 0 -mmq 30 -mbq 0 -orientated -distance $data_matrix_scale -minCTdepth 1 -alignmentType Center ";
			
			$java_cmd .= "-feature $feature_dir"."$feature_file ";
			$java_cmd = $java_cmd."-L $feature_dir"."$region_file ";	
			my $gch_file = $data_dir."$prefix.$genome_assembly.$sample.gch.$data_matrix_scale.txt";
			$java_cmd = $java_cmd."-gchFile $gch_file ";
			my $wcg_file = $data_dir."$prefix.$genome_assembly.$sample.wcg.$data_matrix_scale.txt";
			$java_cmd = $java_cmd."-wcgFile $wcg_file ";
			my $hcg_file = $data_dir."$prefix.$genome_assembly.$sample.hcg.$data_matrix_scale.txt";
			$java_cmd = $java_cmd."-hcgFile $hcg_file ";
			$header .= "$java_cmd\n";
			$prefix .= ".$sample";
			my $R_cmd_plot = "/home/uec-00/shared/production/software/R-2.14.1/R-2.14.1/bin/R --no-restore --no-save --args wd=$plot_dir prefix=$prefix gchfn=$gch_file hcgfn=$hcg_file step=$plot_step scale=$plot_scale axistep=$plot_axis_step smooth=$smooth < $R_plot_script \n";
			$header .= "$R_cmd_plot ";			
			my $outfileName = "submitPBS" . $count;
			open(OUT, ">$outfileName") or die "can not open file:$!";
			print OUT $header;
			#print $header;
			close OUT;
			system("qsub $outfileName");
			$count++;
			print "$feature_file is submitted\n";
			unlink $outfileName;
#			print "$header\n";
		}
		}

		
	}
closedir(DH1);
closedir(DH);
